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ABSTRACT 

Radiative transfer simulations are now at the forefront of numerical astrophysics. They are 
becoming crucial for an increasing number of astrophysical and cosmological problems; at 
the same time their computational cost has come to the reach of currently available computa- 
tional power. Further progress is retarded by the considerable number of different algorithms 
(including various flavours of ray-tracing and moment schemes) developed, which makes the 
selection of the most suitable technique for a given problem a non-trivial task. Assessing the 
validity ranges, accuracy and performances of these schemes is the main aim of this paper, for 
which we have compared 1 1 independent RT codes on 5 test problems: (0) basic physics, (1) 
isothermal H II region expansion and (2) H II region expansion with evolving temperature, 
(3) I-front trapping and shadowing by a dense clump, (4) multiple sources in a cosmological 
density field. The outputs of these tests have been compared and differences analyzed. The 
agreement between the various codes is satisfactory although not perfect. The main source 
of discrepancy appears to reside in the multi-frequency treatment approach, resulting in dif- 
ferent thicknesses of the ionized-neutral transition regions and the temperature structure. The 
present results and tests represent the most complete benchmark available for the develop- 
ment of new codes and improvement of existing ones. To this aim all test inputs and outputs 
are made publicly available in digital form. 

Key words: H II regions — ISM: bubbles — ISM: galaxies: halos — galaxies: high-redshift — 
galaxies: formation — intergalactic medium — cosmology: theory — radiative transfer — meth- 
ods: numerical 



1 INTRODUCTION 

Numerous physical problems require a detailed understanding of 
radiative transfer (RT) of photons in different environments, rang- 



ing from intergalactic and interstellar medium to stellar or plane- 
tary atmospheres. In particular, a number of problems of cosmo- 
logical interest cannot be solved without incorporating RT calcula- 
tions, e.g. modeling and understanding of the Lya forest, absorp- 
tion lines in spectra of high-z quasars, radiative feedback effects, 
the reionization of the intergalactic medium (IGM) and star forma- 
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tion, just to mention a few (see e.g. ICiardi & Ferra 73 l2004 for a 
recent review on some of these topics). In many of these situations, 
the physical conditions are such that the gas in which photons prop- 
agate is optically thick; also, the geometry of the problem often is 
quite complex. As a consequence, approaches relying on optically 
thin or geometrical approximations yield unsatisfactory (and some- 
times incorrect) results. 

The RT equation in 3D space has seven (three spatial, two 
angular, one frequency, one time) dimensions. Although in spe- 
cific cases certain kinds of symmetry or approximations can be 
exploited, leading to a partial simplification, most problems of as- 
trophysical and cosmological interest remain very complex. For 
this reason, although the basic physics involved is well understood, 
the detailed solution of the complete radiative transfer equation is 
presently beyond the available computational capabilities. In addi- 
tion, the technical implementation of the RT equation in numerical 
codes is a very young and immature subject in astrophysics. 

RT approaches have been attempted in the past to study spe- 
cific problems as the light curves of supernovae, radiation from 
proto-stellar and active galactic nuclei accretion disks, line radia- 
tion from collapsing molecular clouds, continuum photon escape 
from galaxies and the effects of dust obscuration in galaxies. Typi- 
cally, either the low dimensionality of these approaches or the sim- 
plified physics allowed a numerical treatment that resulted in a suf- 
ficiently low computational cost given the available machines. As 
this last constraint has become less demanding, the number and 
type of desirable applications has expanded extremely rapidly, par- 
ticularly spreading to galaxy formation and cosmology fields. 

Following this wave of excitement, several groups then at- 
tacked the problem from a variety of perspectives by using com- 
pletely different, independent, and dedicated numerical algorithms. 
It was immediately clear that the validation and assessment of the 
various codes was crucial in order not to waste (always) limited 
computational and human resources. At that point the community 
faced the problem that, in contrast with e.g. gas-dynamic studies, 
very few simple radiative transfer problems admit exact analytical 
solutions that could be used as benchmarks. Just a few years ago 
only a few RT codes were available and these were still mostly 
in the testing/optimization phase, and therefore lacking the neces- 
sary degree of stability required to isolate truly scientific results 
from uncertainties due to internal programming bugs. In the last 
few years the subject has rapidly changed. Not only has the re- 
liability of existing codes matured, a new crop of codes based on 
novel techniques has been developed, making a comparison of them 
timely. The time is now ripe to do this comparison project in a fairly 
complete and meaningful form and this paper presents the detailed 
outcome of our effort. 

The aim of the present comparison is to determine the type 
of problems the codes are (un)able to solve, to understand the 
origin of the differences inevitably found in the results, to stimu- 
late improvements and further developments of the existing codes 
and, finally, to serve as a benchmark to testing future ones. 
We therefore invite interested RT researchers to make use of 
our results, including tests descriptions, input and output data, 
all of which can be found in digital form at the project web- 
site jhttp : / /www . mpa-garching . mpg . de/tsu3/ | At this 
stage our interest is not focused on the performances of the codes 
in terms of speed and optimization. 

This project is a collaboration of most of the community and 
includes a wide range of different methods (e.g. various versions 
of ray-tracing, as well as moment schemes), some of which have 
already been applied to study a variety of astrophysical and cosmo- 



logical problems. The comparison is made among 1 1 independent 
codes, each of which is described concisely in § |2| (and more ex- 
tensively in the corresponding methodology paper, whenever avail- 
able). We would like to emphasize that the interaction among the 
various participating groups has been characterized by a very con- 
structive and scientifically honest attitude, and resulted in many im- 
provements of the codes. 

Here we present the results from a set of tests on fixed den- 
sity fields, both homogeneous and inhomogeneous, which verify 
the radiative transfer methods themselves. In a follow-up paper we 
plan to discuss the direct coupling to gas-dynamics and compare 
the results on a set of several radiative hydrodynamics problems. 



2 THE CODES 

In this section we briefly describe the eleven radiative transfer 
codes which are taking part in this comparison project. The de- 
scriptions point to the more detailed methodology papers whenever 
available. All the codes, their names, authors and current features 
are summarized in Table [I] 

2.1 C^-Ray: Photon-conserving transport of ionizing 
radiation (G. Mellema, I. Iliev, P. Shapiro, M. Alvarez) 

C^-Ray is a grid-based short characteristics ray-tracing code which 
is photon-conserving and causally traces the rays away from the 
ionizing sources up to each cell. Explicit photon-conservation is 
assured by taking a finite-volume approach when calculating the 
photoionization rates, and by using time-averaged optical depths. 
The latter property allows for integration time steps much larger 
than the ionization time scale, which leads to a considerable speed- 
up of the calculation and facilitates the coupling of our code to 
gasdynamic evolution. The code is described and tested in detail in 
Mellema et al. (2006a). 

The frequency dependence of the photoionization rates and 
the photoionization heating rates is dealt with by using frequency- 
integrated rates, stored as functions of the optical depth at the ion- 
ization threshold. In its current version the code includes only hy- 
drogen and does not include the effects of helium, although they 
could be added in a relatively straightforward way. 

The transfer calculation is done using short characteristics, 
where the optical depth is calculated by interpolating values of grid 
cells lying closer to the source. Because of the causal nature of the 
ray-tracing, the calculation cannot easily be parallelized through 
domain decomposition. However, using OpenMP the code is effi- 
ciently parallelized over the sources. The code is currently used for 
large-scale simula t ions of cosmic reioniz ation and its observability 
Jlliev et alJl 2005bVMellema et al."2006b") on grid sizes up to 406^ 
and up to more than 10^ ionizing sources. 

There are ID, 2D and 3D versions of the code available. It was 
developed to be directly coupled with hydrodynamics calculations. 
The large time steps allowed for the radiative transfer enable the 
use of the hydrodynamic time step for evolving the combined sys- 
tem. The first gasdyn amic application of our code is presented in 
iMellema et al] j2005l) . 

2.2 OTVET: Optically-Thin Variable Eddington Tensor 
Code (N. Gnedin, T. Abel) 

The Op tically Thin Variable E ddington Tensor (OTVET) approxi- 
mation iOnedin & Abell2001l) is based on the moment formulation 
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Table 1. Participating codes and their current features. 



Code (Authors) Grid Parallehzation gasdynamics Hehum Rec. radiation 



C^-Ray (G. Mellema, 1. Ihev, P. Shapiro, M. Alvarez) 


fixed/AMR 


shared 


yes 


no 


no 


OTVET (N. Gnedin, T. Abel) 


fixed 


shared 


yes 


yes 


yes 


CRASH (A. MaseUi, A. Ferrara, B. Ciardi) 


fixed 


no 


no 


yes 


yes 


RSPH (H. Susa, M. Umemura) 


no grid, particle-based 


distributed 


yes 


no 


no 


ART (T. Nakamoto, H. Susa, K. Hiroi, M. Umemura) 


fixed 


distributed 


no 


no 


yes 


FTTE (A. Razoumov) 


fixed/AMR 


no 


yes 


yes 


yes 


Simplex (J. Ritzerveld, V. Icke, E.-J. Rijkhorst) 


unstructured 


no 


no 


no 


yes 


Zeus-MP (D. Whalen, M. Norman) 


fixed 


distributed 


yes 


no 


no 


FLASH-HC (E.-J. Rijkhorst, T. Plewa, A. Dubey, G. Mellema) 


fixed/AMR 


distributed 


yes 


no 


no 


IFT (M. Alvarez, R Shapiro) 


fixed/AMR 


no 


no 


no 


no 


Coral (I. IHev, A. Raga, G. Mellema, P. Shapiro) 


AMR 


no 


yes 


yes 


no 



of the radiative transfer equation: 



c dt dx^ 
o dF^ d 



(1) 



where Ei^ and F^ are the energy density and the flux of radiation 
respectively, is the absorption coefficient, is the source func- 
tion, and hl^ is a unit trace tensor normally called the Eddington 
tensor. 

It is important to underscore that the source function Su{x) 
is considered to be an arbitrary function of position, so that it may 
contain both the delta-function contributions from any number of 
individual point sources and the smoothly varying contributions 
from the diffuse sources. 

Equations (ij form an open system of two partial differen- 
tial equations, because the Eddington tensor can not be determined 
from them. In the OTVET approximation the Eddington tensor is 
computed from all sources of radiation as if they were optically 
thin, hi^ = Pi^/Ti Pi\ where 



Ax' - x\){x^ - x{) 



{x-xiY 



where is the mass density of the sources e.g. stars. 

Thus, the OTVET approximation conserves the number den- 
sity of photons (in the absence of absorption) and the flux, but may 
introduce an error in the direction of the flux propagation. It can 
be shown rigorously that the OTVET approximation is exact for a 
single point-like source and for uniformly distributed sources, but 
is not exact in other cases. 

While it is difficult to prove rigorously, it appears that the 
largest error is introduced for the case of two sources, one much 
stronger than the other. In that case the H II region around the strong 
source is modeled highly precisely, but the shape of the H II region 
around the weaker source (before the two H II regions merge) be- 
comes ellipsoidal with the deviation from the spherical symmetry 
never exceeding 17% (1/6) in any direction. 



2.3 CRASH: Cosmological RAdiative transfer Scheme for 
Hydrodynamics (A. Maselli, A. Ferrara, B. Ciardi) 

CRASH is a 3D ray-tracing radiative-transfer code based on the 
Monte Carlo (MC) technique for sampling distribution functions. 



The grid-based algorithm follows the propagation of the ionizing 
radiation through an arbitrary H/He static density field and calcu- 
lates the time-evolving temperature and ionization structure of the 
gas. 

The MC approach to RT requires that the radiation field is dis- 
cretized into photon packets. The radiation field is thus reproduced 
by emitting packets, according to the configuration under analy- 
sis, and by following their propagation accounting for the opacity 
of the gas. For each emitted photon packet, the emission location, 
frequency and propagation direction are determined by randomly 
sampling the appropriate probability distribution functions (PDFs), 
which are assigned as initial conditions. It is possible to include 
an arbitrary number of point/extended sources and/or diffuse back- 
ground radiation in a single simulation. This approach thus allows 
a straightforward and self-consistent treatment of the diffuse radia- 
tion produced by H/He recombinations in the ionized gas. 

The relevant radiation-matter interactions are accounted for 
during the photon packet's propagation. At each cell crossed, each 
packet deposits a fraction of its photon content according to the 
cell's opacity which determines the absorption probability. Once 
the number of photons absorbed in the cell is calculated, we find 
the effect on the temperature and on the ionization state of the gas, 
by solving the discretized non-equilibrium chemistry and energy 
equations. Recombinations, collisional ionizations and cooling are 
treated as continuous processes. 

The detailed description of the CRASH imple mentation is 
given in Ciardi et alJ <20Cnh and lMaselli etalJ J2003 V and an im- 
proved algorithm for de aling with a background d iffuse ionizing 
radiation is described in iMaselli & Ferrara i2005h . The tests de- 
scribed in this paper have been performed using the most updated 
version of the code (Maselli et al. 2003). 

The code has been primarily developed to study a number of 
cosmological problems, such as hydrogen and helium reionization, 
the physical state of the Lya forest, the escape fraction of Lyman 
continuum photons from galaxies, the diffuse Lya emission from 
recombining gas. However, its flexibility allows applications that 
could be relevant to a wide range of astrophysical problems. The 
code architecture is sufficiently simple that additional physics can 
be easily added using the algorithms already implemented. For ex- 
ample, dust absorption/re-emission can be included with minimum 
effort; molecular opacity and line emission, although more compli- 
cated, do not represent a particular challenge given the numerical 
scheme adopted. Obviously, were such processes added, the com- 
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putational time could become so long that parallelization would 
be necessary. This would be required also when CRASH will be 
coupled to a hydrodynamical code to study the feedback of photo- 
processes onto the (thermo-)dynamics of the system. 



2.4 RSPH: SPH coupled with radiative transfer (H. Susa, M. 
Umemura) 

The Radiation-SPH scheme is designed to investigate the formation 
and evolution of the first generation objects at z ^ 10 (Susa 2006), 
where the radiative feedback from various sources play important 
roles. The code can compute the fraction of chemical species e, 
H+, H, H~, H2, and by fully implicit time integration. It also 
can deal with multiple sources of ionizing radiation, as well as the 
radiation in the Lyman- Werner band. 

Hydrodynamics is calculated by the Smoothed Parti cle Hydro- 
dynam ics (SPH) method. We use the version of SPH by'Umemura' 
Jl993h with the modification by Steinmetz & Mueller ( 1993), and 
we also adopt the particle resizing formalism by Thacker et al. 
(Eobo,) . In the present version, we do not use the entropy formal- 
ism. 

The non-equilibrium chemistry and radiative cooling for 
primordial gas are calculated by the code developed by 
[Susa & Kitava ma ( 200 Q), where H2 co oling and reaction rates are 
mostly taken from Galli & Pallj Jl998h . 

As for the photoionization process, we employ the on-the-spot 
approximation (S Ditzeiiil978 ). We solve the transfer of ionizing 
photons directly from the source but do not solve the transfer of dif- 
fuse photons. Instead, it is assumed that the recombination photons 
are absorbed in the neighbourhood of the spatial position where 
they are emitted. The absence of source terms in this approxima- 
tion greatly simplifies the radiation transfer equation. Solving the 
transfer equation reduces to the determination of the optical depth 
from the source to every SPH particle. 

The optical depth is integrated utilizing the neighbour 
lists of SPH par t icles. It is similar to the code described in 
ISusa & Umemur j fe004). but now we can deal with multiple point 
sources. In our new scheme, we do n ot create so many grid p oints 
on the light ray as the previous code iSusa & Umemural2004) did. 
Instead, we just create one grid point per SPH particle in its neigh- 
bor. We find the 'upstream' particle for each SPH particle on its 
line of sight to the source. Then the optical depth from the source 
to the SPH particle is obtained by summing up the optical depth at 
the 'upstream' particle and the differential optical depth between 
the two particles. 

The code is already parallelized using the MPI library. The 
computational domain is divided by the Orthogonal Recursive Bi- 
section method. The parallelization method for radiation trans- 
fer p^H^_simil^to the Multipl e Wave Front m ethod developed 
by 'Nakamoto et al.' ("2001) and iHeinemann e t al. (2005), but it is 
changed to fit the RSPH code. The details are described in Susa 
The code also is able to handle gravity with a Barnes-Hut 
tree, which is also parallelized. 



2.5 ART: Authentic Radiative Transfer with Discretized 
Long Beams (T. Nakamoto, H. Susa, K. Hiroi, M. 
Umemura) 

ART is a grid based code designed to solve the transfer of radia- 
tion from point sources as well as diffuse radiation. On the photo- 
ionization problem, ART solves time-dependent ionization states 



and energies. Hydrodynamics is not incorporated into the current 
version. 

In the first step of the ART scheme, a ray, on which photons 
propagate, is cast from the origin by specifying the propagation an- 
gles. When the distance from the ray to a grid point is smaller than 
a certain value, a segment is located at the grid point. A collection 
of all the segments along the ray is considered to be a decomposi- 
tion of the ray into segments. The radiative transfer calculation is 
sequentially done from the origin towards the downstream side on 
each segment. From one segment to another, the optical depth is 
calculated and added. Finally, changing the angle and shifting the 
origin, we can obtain intensities at all the grid points directed along 
all the angles. 

ART has two versions for the integration of intensities over 
the angle. The first is a simple summation of intensities with fi- 
nite solid angles. This scheme is fit to diffuse radiation cases. The 
second one, designed to fit point sources, uses a solid angle of the 
source at the grid point to evaluate the dilution factor. Multiplying 
the representing intensity at the grid point and the dilution factor, 
we can obtain the integration of the intensity over the angle. Gen- 
erally, the radiation field can be divided into two parts: one is the 
direct incident radiation from point sources and the other is the dif- 
fuse radiation. ART can treat both radiation fields appropriately by 
using two schemes simultaneously. 

The integration over the frequency is done using the one- 
frequency method, which is similar to the six-frequency method 
devised by Nakamoto et al. (2001). Since in our present problems 
both the spectrum of the source and the frequency-dependence of 
the absorption coefficient are known in advance, once we calcu- 
late the optical depth at the Lyman limit frequency, the amount of 
absorption at any frequency can be obtained without carrying out 
integration along the ray for the different frequency. 

The parallelization of ART can be done based not only on 
the angle-frequency decomposition but also on the spatial domain 
decomposition using the Multiple Wave Front (MWF) method 
iNakamoto et all200lh. 



2.6 FTTE: Fully threaded transport engine (A. Razoumov) 

FTTE is a new method for transport of both diffuse and point 
source radiation on refined grids developed at Oak Ridge National 
Laboratory. The diffuse part of the solver has been described in 
Razoumov & Cardall ( 2005). Transfer around point sources is done 
in a separate module acting on the same fully threaded data struc- 
ture (3D fields of density, temperature, etc.) as the diffuse part. The 
point sour ce algorithm is an extens ion of the adaptive ray- splitting 
scheme of I Abel & WandelJ i2002h to a model with variable grid 
resolution, with all discretization done on the grid. Sources of ra- 
diation can be hosted by cells of any level of refinement, although 
usually in cosmological applications sources reside in the deepest 
level of refinement. Around each point source we build a system 
of radial rays which split either when we move further away from 
the source, or when we enter a refined cell, to match the local min- 
imum required angular resolution. Once any radial ray is refined 
it stays refined (even if we leave the high spatial resolution patch) 
until further angular refinement is necessary. 

All ray segments are stored as elements of their host cells, and 
actual transport just follows these interconnected data structures. 
Whenever possible, an attempt is made to have the most generic 
ray pattern possible, as each ray pattern needs to be computed only 
once. The simplest example is an unigrid calculation (no refine- 
ment), where there is a single ray pattern which can be used for all 
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sources. Another computationally trivial example is a collection of 
halos within refined patches with the same local grid geometry as 
seen from each source. 

For multiple sources located in the same H II region, we can 
merge their respective ray trees when the distance from the sources 
to a ray segment far exceeds the source separatio n, and the opt ical 
depth to both sources is negligible - see RazoumovetaDIZQOa) for 
details. 

The transport quantity in the diffuse solver is the specific (per 
unit frequency) intensity, whereas in the point source solver it is 
the specific photon luminosity - the number of photons per unit 
frequency entering a particular ray segment per unit time. For dis- 
cretization in angles both modules use the HEALPix algorithm 
iGorski et all2002 ). dividing the entire sphere into 12 x 4"^"^ equal 
area pixels, where n — 1, 2, ... is the local angular resolution. For 
point source transfer, n is a function of the local grid resolution, 
the physical distance to the source, and the local ray pattern in a 
cell. As we go from one ray segment to another, in each cell we 
accumulate the mean diffuse intensity and all photo-reaction rates 
due to point source radiation. Fo r rate equations, we u se the time- 
dependent chemistry solver from'Anninos e t alJil997h . 

Currently there is only a serial version of the code, although 
there is a project underway at the San Diego Supercomputing Cen- 
ter to parallelize the diffuse part of the algorithm angle-by- angle. 
Even in the serial mode the code is very fast, limited in practice by 
the memory available to hold ray patterns. We have run various test 
problems up to grid sizes 256^ with five levels of refinement (8192^ 
effective spatial resolution), and up to the angular resolution level 
n = 13. 

2.7 Simplex : Radiative Transfer on Unstructured Grids (J. 
Ritzerveld, V. Icke, E.-J. Rijkhorst) 

Simplex jRitzerveld et al."2003') is a mesoscopic particle method, 
using unstructured Lagrangian grids to solve the Boltzmann equa- 
tion for a photon gas. It has many similarities with Lattice Boltz- 
mann Solvers, which are used in complex fluid flow simulations, 
with the exception that it uses an adaptive grid based on the cri- 
terion that the local grid step is chosen to correlate with the mean 
free path of the photon. In the optimal case, this last modification 
results in an operation count for our method which does not scale 
with the number of sources. 

More specifically, SimpleX does not use a grid in the usual 
sense, but has as a basis a point distribution, which follows the 
medium density, or the opacity. Given a regular grid with medium 
density values, we use a Monte Carlo method to sample our points 
according to this density distribution. We construct an unstructured 
grid from this point distribution by using the Delaunay tessellation 
technique. This recipe was intentionally chosen this way to ensure 
that the local optical mean free paths correlate linearly with the line 
lengths between points. This way, the radiation-matter interactions, 
which determine the collision term on the rhs of the Boltzmann 
equation., can be incorporated by introducing a set of 'interaction 
coefficients' {c^}, one for each interaction. These coefficients are 
exactly the linear correlation coefficients between the optical mean 
free path and the local line lengths. 

The transport of radiation through a medium can subsequently 
be defined and implemented as a walk on this resultant graph, with 
the interaction taking place at each node. The operation count of 
this resultant method is O {N^^^^^^, in which N is the number of 
points, or resolution, and m is the dimension. This is ^dependent 
of the number of sources, which makes it ideal to do large scale 



reionization calculations, in which a large number of sources is 
needed. 

The generality of the method's setup defines its versatility. 
Boltzmann-like transport equations describe not only the flow of 
a photon gas, but also that of a fluid or of a plasma. It is therefore 
straightforward to define a SimpleX method which solves both the 
radiative transfer equations and the hydrodynamics equations self- 
consistently. At the same time, we are in the process of making 
a dynam ic coupling of Sini pleX with the gri d based hydro codes 
FLASH JFrvxell etall200Ch and GADGET-2 iSDringei20o3) . The 
code easily runs on a single desktop machine, but will be paral- 
lelized in order to accommodate this coupling. 

For this comparison project, we used a SimpleX method set 
up to do cosmological radiative transfer calculations. The result 
is a method which is designed to be photon conserving, updat- 
ing the ionization fraction and the resultant exact local optical 
depth dynamically throughout the simulation. Moreover, diffuse 
recombination radiation, shown to be quite important for the over- 
all result fe.g. iRitzerveldEoO^ . can readily be implemented self- 
consistently and without loss of computational speed. 

For now, we do not solve the energy equation, and do simula- 
tions H-only. We ignore the effect of spectral hardening, by which 
we can analytically derive the fraction of the flux above the Lyman 
limit, given that the source radiates as a black body. We account 
for H-absorption by using a blackbody averaged absorption coeffi- 
cient. The final results are interpolated from the unstructured grid 
cells onto the 128^ datacube to accommodate this comparison. 

2.8 ZEUS-MP with radiative transfer (D. Whalen, M. 
Norman) 

The ZEUS-MP hydrocode in this comparison project solves ex- 
plicit finite-difference approximations to Euler's equations of fluid 
dynamics together with a 9- species reactive network that utilizes 
photoionization rate coefficients computed by a ray-casting radia- 
tive transfer module dWhalen & Nor man 2006). Ionization fronts 
thus arise as an emergent feature of reactive flow and radiative 
transfer in our simulations and are not tracked by computing equi- 
libria positions along lines of sight. The hydrodynamical vari- 
ables (p, e, and the pVi) are updated term by term in operator- 
split and directionally- split substeps, with a given substep incor- 
porating the partial update from the previous substep. The gradient 
(force) terms in the Euler equations are computed in source rou- 
tines and th e divergence terms are calculated in advection routines 
JStone & Nor man 1992). 

The primordial species added to ZEUS-MP (H, H+, He, 
He^, He^^, H~, HJ, H2, and e~) are evolved by nine addi- 
ti onal continuity equa tions and the nonequilibrium rate equations 
of lAnninos et alJ il99 7). The divergence term for each species is 
evaluated in the advection routines, while the other terms form a 
reaction network that is solved separately from the source and ad- 
vective updates. Although the calculations performed for this com- 
parison study take the gas to be hydrogen only, in general we se- 
quentially advance each n^ in the network, building the i*^ species' 
update from the i - 1 (and earlier) updated species while apply- 
ing rate coefficients evaluated at the current problem time. Charge 
and baryon conservation are enforced at the end of each hydrody- 
namic cycle and microphysical cooling and heating processes are 
included by an isochoric operator- split update to the energy density 
computed each time the reaction network is advanced. The radia- 
tive transfer module computes the photoionization rate coefficients 
required by the reaction network by solving the static equation of 
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transfer recast into flux form along radial rays outward from a point 
source centered in a spherical-coordinate geometry. The number of 
ionizations in a zone is calculated in a photon-conserving manner 
to be the number of photons entering the zone minus the number 
exiting. 

The order of execution of the algorithm is as follows: first, the 
radiative transfer module is called to calculate ionization rates in 
order to determine the smallest heating/cooling time on the grid. 
The grid minimum of the Courant time is then computed and the 
hydrodynamics equations are evolved over the smaller of the two 
timescales. Next, the shortest chemistry timescale of the grid is cal- 
culated 

Atchem = 0.1 — (2) 
Tie 

which is formulated to ensure that the fastest reaction operating 
at any place or time in the problem governs the maximum time 
by which the reaction network may be accurately advanced. The 
species concentrations and gas energy are then advanced over this 
timestep, the transfer module is called again to compute a new 
chemistry timestep, and the network and energy updates are per- 
formed again. The n^ and energy are subcycled over successive 
chemistry timesteps until the hydrodynamical timestep has been 
covered, at which point full updates of velocities, energies, and 
densities by the source and advection routines are computed. A new 
hydrodynamical timestep is then determined and the cycle repeats. 

This algorithm has been extensively tested with a comprehen- 
sive suite of quantitative static and hydrodynamical benchmarks 
complementing those appearing in this paper. The tests are de- 
scribed in detail in lWhalen & NormanI (l2006h 



2.9 FLASH-HC: Hybrid Characteristics (E.-J. Rijkhorst, T. 
Plewa, A. Dubey, G. Mellema) 

The Hybrid Chara cteristics (HC) method feiikhorsJ l2005l: 
iRiikhorst et all2006h is a three-dimensional radiative transfer algo- 
rithm designed specifically for use with parallel adaptive mesh re- 
finement (AMR) hydrodynamics codes. It introduces a novel form 
of ray tracing that can neither be classified as long, nor as short 
characteristics. It however does apply the underlying principles, 
i.e. efficient execution through interpolation and parallelizability, 
of both these approaches. 

Primary applications of the HC method are radiation hydrody- 
namics problems that take into account the effects of photoioniza- 
tion and heating due to point sources of radiation. T he method is 
imp lemented into the hydrodynamics package Flash iFrvxell et alJ 
|2000). The ionization, heating, and cooling processes are modeled 
using the Doric package (Frank & Mellema 1994). Upon compari- 
son with the long characteristics method, it was found that the HC 
method calculates shadows with similarly high accuracy. Although 
the method is developed for problems involving photoionization 
due to point sources, the algorithm can easily be adapted to the 
case of more general radiation fields. 

The Hybrid Characteristics algorithm can be summarized as 
follows. Consider an AMR hierarchy of grids that is distributed 
over a number of processors. Rays are traced over these different 
grids and must, to make this a parallel algorithm, be split up into 
independent ray sections. Naturally these sections are in the first 
place defined by the boundaries of each processor's sub-domain, 
and in the second place by the boundaries of the grids contained 
within that sub-domain. At the start of a time step, each proces- 
sor checks if its sub-domain contains the source. The processor 



that owns the source stores its grid and processor id and makes it 
available to all other processors. Then each processor ray traces the 
grids it owns to obtain local column density contributions. Since in 
general rays traverse more than one processor domain, these local 
contributions are made available on all processors through a global 
communication operation. By interpolating and accumulating all 
contributions for all rays (using a so called 'grid-mapping') the to- 
tal column density for each cell is obtained. The coefficients used in 
the interpolation are chosen such that the exact solution for the col- 
umn density is retrieved when there are no gradients in the density 
distribution. Tests with a 1/r^ density distribution resulted in er- 
rors < 0.5% in the value for the total column density as compared 
to a long characteristics method. 

Once the column density from the source up to each cell face is 
known, the ionization fractions and t emperature can be computed . 
For this we use the Doric package fsee lMellema & Lu ndavist 2 00^ 
Frank & Mellema 1994). These routines calculate the photo- and 
collisional ionization, the photoheating, and the radiative cooling 
rate. Using an analytical solution to the rate equation for the ion- 
ization fractions, the temperature and ionization fractions are found 
through an iterative process. Since evaluating the integrals for the 
photoionization and heating rate is too time consuming to perform 
for every value of the optical depth, they are stored in look-up ta- 
bles and are interpolated when needed. 

The hydrodynamics and ionization calculations are coupled 
through operator splitting. To avoid having to take time steps 
that are the minimum of the hydrodynamics, ionization, and heat- 
ing/cooling time scales, the fact that the equations for the ioniza- 
tion and heating/cooling can be iterated to convergence is used. 
This means that the only restriction on the time step comes from 
the hydrodynamics (i.e. the Courant condition). Note however that, 
when one needs to follow R-type ionization fronts, an additional 
time step constraint is used to find the correct propagation velocity 
for this type of front. 

An assessm ent of the parallel performance of the HC method 
was presented by RiikhorsJ i2005h . It was found that the ray trac- 
ing part takes less time to execute than other parts of the calcu- 
lation (e.g. hydrodynamics and AMR.). Tests involving randomly 
distributed sources show that the algorithm scales linearly with the 
number of sources. Weak scaling tests, where the amount of work 
per processor is kept constant, as well as strong scaling tests, where 
the total amount of work is kept constant, were performed as well. 
By carefully choosing the amount of work per processor it was 
shown that the HC algorithm scales well for at least '^100 pro- 
cessors on an SGI Altix, and ^ 1000 processors on an IBM Blue- 
Gene/L system. 

The HC method will be made publicly available in a future 
Flash release. 



2.10 IFT: Ionization Front Tracking (M. Alvarez and P. 
Shapiro) 

We have developed a ray-tracing code which explicitly follows the 
progress of an ionization front (I-front) around a point source of 
ionizing radiation in an arbitrary three-dimensional density field 
of atomic hydrogen. Because we do not solve the non-equilibrium 
chemical and energy rate equations and use a simplified treatment 
of radiative transfer, our method is extremely fast. While our code is 
currently capable of handling only one source, we are in the process 
of generalizing it to handle an arbitrary number of point sources. 
More detailed descriptions of various aspects of this code can be 
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found in § 5.1.3 of lMellema etalJJ2006j) and §3 of I Alvarez etall 
J2005h . 

The fundamental assumption we make is that the I-front is 
sharp. Behind the front, the ionized fraction and temperature are 
assumed to take their equilibrium values, while ahead of the front 
they are assigned their initial values. The assumption of equilib- 
rium behind the I-front is justified because the equilibration time 
on the ionized side of the I-front is shorter than the recombina- 
tion time by a factor of the neutral fraction, which is small behind 
the I-front. Because the progress of the I-front will be different in 
different directions, it is necessary to solve for its time-dependent 
position along rays that emanates from the source, the angular ori- 
entation of wh ich are chosen to lie at the centers of HEALPixels^ 
iGorski et al 12005 ). Typically, we choose a sufficient resolution of 
rays in the sky such that there is approximately one ray per edge 
cell. Along a given ray, we solve the fully-relativistic equation for 
the propagation of the I-front ( Sh apiro et al 2005.) : 

dR ^ cQ{R) 

dt Q{R) +47rR^cn{R)' ^ ^ 

where Q{R) is the ionizing photon luminosity at the surface of 
the front, R is the distance along the ray, and c is the speed of 
light. This equation correctly accounts for the finite travel time of 
ionizing photons, i.e. as Q{R) — > oo, dR/dt — > c. The arrival 
rate of ionizing photons is given by 

nR 

Q(R) = - 47r / aB{r)n{rydr, (4) 
Jo 

where Q* is the ionizing photon luminosity of the source, n(r) is 
the density along the ray, and asir) is the "case B" recombination 
coefficient along the ray. The value of as (r) is determined by the 
equilibrium temperature, which varies along the ray. 

A brief outline of the algorithm is as follows. First, we inter- 
polate the density from the grid to discrete points along each ray, 
where the spacing between points along each ray is approximately 
the same as the cell size of the grid. Next, we compute the equi- 
librium profile of ionized fraction and temperature along each ray 
by moving outward from the source, using the equilibrium neutral 
hydrogen density of the previous zones to attenuate the flux to each 
successive zone. Equation (|3} is then solved for each ray, which 
gives the I-front position in that direction. The values of ionized 
fraction and temperature along each ray are set to their equilibrium 
values inside of the I-front, and set to their initial values on the out- 
side. Finally, the ionized fraction and temperature are interpolated 
back from the rays to the grid. 



2.11 CORAL (I. Iliev, A. Raga, G. Mellema, P. Shapiro) 

CORAL is a 2-D, axisymmetric Eulerian fluid dynamics adap- 
tive mesh refinem ent (AMR) code (see iMellema et alJ Il998li 
IShapiro et all2004 and references therein for detailed description). 
It solves the Euler equations in their conservative finite-volume 
form using the second-order method of van Leer flux- splitting, 
which allows for correct and precise treatment of shocks. The grid 
refinement and de-refinement criteria are based on the gradients of 
all code variables. When the gradient of any variable is larger than 
a pre-defined value the cell is refined, while when the criterion for 
refinement is not met the cell is de-refined. 



^ http://healpix.jpl.nasa.gov 
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Figure 1. Legend for the line plots. 

The code follows, by a semi-implicit method, the non- 
equilibrium chemistry of multiple species (H, He, C II- VI, N I- 
yi, O I- VI, Ne I- VI, and S II- VI) and the corresponding cooling 
jRaga et al. 1997; Mellema et al. 1998), as well as Compton cool- 
ing. The photoheating rate used is a sum of the photoionization 
heating rates for H I, He I and He II. For computational efficiency 
all heating and cooling rates are pre-computed and stored in tables. 
The microphysical processes - chemical reactions, radiative pro- 
cesses, transfer of radiation, heating and cooling - are implemented 
though the standard approach of operator- splitting (i.e. solved each 
time- step, side-by- side with the hydrodynamics and coupled to it 
through the energy equation). The latest versions of the code also 
include the effects of an external gravity force. 

Currently the code uses a black-body, or a power-law ionizing 
source spectra, although any other spectrum can be accommodated. 
Radiative transfer of the ionizing photons is treated explicitly by 
taking into account the bound-free opacity of H and He in the pho- 
toionization and photoheating rates. The photoionization and pho- 
toheating rates of H I, He I and He II are pre-computed for the given 
spectrum and stored in tables vs. the optical depths at the ionizing 
thresholds of these species, which are then used to obtain the total 
optical depths. The code correctly tracks both fast (by evolving on 
an ionization timestep. At ^ uh /uh) and slow I-fronts. 

The code has been tested extensively and has been applied 
to many astrophysical problems, e.g. photo evaporation of clumps 
in planetary nebulae iMellema et alJ ll998D . cosmological mini- 
halo photoevapo ration during reionization fshauiro et al. 2 0041 : 
llliev et alJ20053) . and studies of the radiative feedback from propa- 
gating ionization fronts on dense clumps in Damped Lyman-a sys- 
tems ( Iliev et al. 2005a) . 



3 TESTS AND RESULTS 

In this section we describe the tests we have performed, along with 
their detailed parameters, geometry and setup, and the results we 
obtained. When constructing these tests, we aimed for the simplest 
and cleanest, but nonetheless cosmologically-interesting problems. 
We designed them in a way which allows us to test and compare 
all the important aspects of any radiative-transfer code. These in- 
clude correct tracking of both slow and fast I-fronts, in homoge- 
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Figure 2. Test 0, part 1: Hydrogen rates, cooling and cross-sections used by the participating codes. 



neous and inhomogeneous density fields, formation of shadows, 
spectrum hardening, and solving for the gas temperature state. In 
a companion paper (Paper II) we will present the tests which in- 
clude the interaction with fluid flows and radiative feedback on the 
gas. For simplicity, in all tests the gas is assumed to be hydrogen- 
only. The gas density distribution is fixed. Finally, in order to be 
included in the comparison, each test had to be done by at least 
three or more of the participating codes. 

Figure Q provides a legend allowing the reader to identify 
which line corresponds to which code in the figures throughout the 
paper. The images we present are identified in the corresponding 
figure caption. 

All test problems are solved in three dimensions (3D), with 
grid dimensions 128^ cells. For the ID (Zeus-MP) and 2D (Coral) 
codes the data is interpolated on the same- size 3D grid and ana- 
lyzed the same way. Unless otherwise noted, the sources of ionizing 
radiation are assumed to have a black-body spectrum with effective 
temperature Tgff = 100, 000 K, except for Test 1 where the as- 
sumed spectrum is monochromatic with all photons having energy 
hu — 13.6 eV, the ionization threshold of hydrogen. In all tests the 
temperature is allowed to vary due to atomic heating and cooling 
processes in accordance with the energy equation, again with the 
exception of Test 1 where we fix the temperature at T = 10^ K. 
The few codes that do not yet include an energy equation use a 
constant temperature value. 



3.1 Test 0: The basic physics 

The solution of the radiative transfer equation is intimately related 
to the ionization and thermal states of the gas. These depend on 
the atomic physics reaction rates, photoionization cross- sections, 
as well as cooling and heating rates used. As there is a variety 
of rates available in the literature, for the sake of clarity we have 
summarized those used in our codes in Tableland plotted them in 
Figure|2| The table columns indicate, from left to right, the name of 
the code and the reference for: case A recombination rate of H II, 
He II and He III; case B recombination rate of H II, He II and He III; 
dielectronic recombination rate of He II; collisional ionization rate 
of H I, He I and He III; collisional ionization cooling rate of H I, 
He I and He II; case A recombination cooling rate of H II, He II and 
He III; case B recombination cooling rate of H II, He II and He III; 
dielectronic recombination cooling rate of He II; collisional exci- 
tation cooling rate of H I, He I and He II; Bremsstrahlung cooling 
rate; Compton cooling rate; cross-section of H I, He I and He II. 
Note that for those codes in which the treatment of He is not in- 
cluded, only the references to the H rates are given. 

The first thing to note in Figure |2l is that although the rates 
come from a wide variety of sources, they largely agree. The main 
differences are in our recombination cooling rates, particularly at 
very high temperatures, beyond the typical range of gas tempera- 
tures achieved by photoionization heating (T < 10^ K). It should 
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Table 2. Rates adopted by the different radiative transfer codes. The columns are, from left to right, name of the code and reference for: case A recombination 
rate (RRA) of H II, He II and He III; case B recombination rate (RRB) of H II, He II and He III; dielectronic recombination rate (DRR) of He II; coUisional 
ionization rate (CIR) of H I, He I and He III; coUisional ionization cooling rate (CICR) of H I, He I and He II; case A recombination cooling rate (RCRA) 
of H II, He II and He III; case B recombination cooling rate (RCRB) of H II, He II and He III; dielectronic recombination cooling rate (DRCR) of He II; 
coUisional excitation cooling rate (CECR) of H I, He I and He II; Bremsstrahlung cooling rate (BCR); Compton cooling rate (CCR); cross-section (CS) of H I, 
He I and He 11. Units of RRA, RRB, DRR and CIR are [cm^s-^, units of CICR, RCRA, RCRB, DRCR, CECR, BCR, CCR are [erg cm^s-^ and units of 
CS are [cm^]. Note that for those codes in which the treatment of He is not included, only the references to the H rates are given. 
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(1997); 

Notes: the codes C^-Ray, Flash-HC and Coral aU share the same nonequilibrium chemistry module (DORIC, developed by G. MeUema), so they have the 
same hydrogen chemistry and heating rates and cross-section, and the same chemistry solver. However, Coral also includes the chemistry of helium and a 

number of metals. 



be noted, however, that e.g. shock-heated gas can reach much 
higher temperatures, in which case the differences in our rates be- 
come very large and caution should be excersised in choosing the 
appropriate rates. However, even at the typical photoionization tem- 
peratures there are differences between the rates by up to factors of 
^ 2. The origin of these discrepancies is currently unclear, and it 
is also unclear which fit to the experimental data is more precise. 
There are also a few cases in which particular cooling rates (e.g. 
Zeus Case B recombination cooling, ART line cooling) are notably 
different from the rest. 

As a next step, we did several numerical experiments to assess 
the offsets between our results that can arise solely based on our 
different chemistry and cooling rates. We did this by implementing 
the rates from several codes representative of the full range of rates 
present above into a single code (1-D version of ART code) and 
running the same problem (Test 2 described below, which is expan- 
sion of an H II region in uniform density gas, see ^ 13. 3 I f or detailed 
test definition and solution features). In all cases the same photoion- 
ization cross-section (the one of ART) is used and only the rates 
are varied. The code OTVET here stands also for the codes C^- 
Ray, CRASH, Flash-HC and Coral, since all these codes have either 
identical, or closely-matching rates. Results are shown in Figures|3] 
and0 In Figure|3]we show the I-front position, n, and the I-front 
velocity, vi (Both quantities are normalized to the analytical solu- 
tions of that problem obtained at temperature T = 10^ K, while the 
time is in units of the recombination time at the same temperature). 
Note that the results for Zeus use the recombination cooling rate of 
ART since its currently-implemented rate is overly- simplified. We 
show the results using the rates of OTVET, FTTE, RSPH, Zeus-MP 
and ART. 

Initially, when the I-front is fast and still far away from reach- 
ing its Stromgren sphere the results for all codes agree fairly well, 
as expected since recombinations are still unimportant. Once re- 




Figure 3. Test 0, part 2: I-front expansion in uniform-density field with 
temperature evolution (same as Test 2 below). Plotted are the I-front po- 
sition (top) and velocity (bottom) all derived from the same code (ID, 
spherically-symmetric version of ART code) and using the photoionization 
cross-section from that code, but with the rest of the microphysics (chem- 
istry and cooling rates) taken from several of the participating codes, as 
indicated by line-type and color. All results are normalized to the analytical 
ones (which assume fixed temperature, T = 10^ K) given in equations J6| 
below (see text for details). 
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Figure 5. Test 0, part 3: Single-zone ionizing up and cooling and recombin- 
ing. 



combinations do become important, at t > tree, the results start di- 
verging. The results for OTVET, FTTE and RSPH remain in close 
agreement, within a fraction of once per cent in the I-front radius 
and within ~ 2% in the I-front velocity. The results using the Zeus 
and ART rates however depart noticeably from the others, by up to 
4% and 6, respectively in radius. The corresponding velocities are 
different even more, by up to factor of 2 at the end, when the 
I-front is close to stationary. In Figure |4] we show the radial pro- 
files of the ionized fraction x, and temperature T, normalized to 
the results using OTVET rates (left), and in absolute units (right). 
The ionized fraction profiles for ART and Zeus again are fairly dif- 
ferent from the rest, by 20-40%, while the rest of the codes agree 
between themselves much better, to less than 10%. Agreement is 
slightly worse close to the ionizing source. In terms of temperature 
profiles the codes agree to better than 10%, with the exception of 
ART, in which case the resulting temperature is noticeably lower, 
by - 20%. 

The reasons for these discrepancies lie largely in differences 
in the recombination rates and the recombination and line cooling 
rates. The line cooling rate of ART is larger than the ones for most 
of the other codes by a factor of 2 - 5 in the temperature range 
10^ — 10^ K, while the recombination rate used by that code is 
about 10-15 percent larger in the same temperature range, thus the 
resulting gas temperature is correspondingly lower. This lower tem- 
perature results in a higher recombination rate and hence the slower 
propagation of the I-front we observed. The reason for the discrep- 
ancy with Zeus is mostly due to its higher recombination rate. This 
again results in a somewhat slower I-front propagation, but not in 
significant temperature differences. Accordingly, for both of these 
codes the neutral gas fraction is significantly higher at all radii. The 
slightly higher recombination cooling of RSPH results in slightly 
lower temperature and proportionally higher neutral fraction, al- 
though both are off by only a few percent, and up to 10% close to 
the ionizing source. 

The last important element in this basic physics comparison 
is to assess the accuracy and robustness of the methods we use for 
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Figure 6. Test 1 (H II region expansion in an uniform gas at fixed temperature): Images of the H I fraction, cut through the simulation volume at coordinate 
z = at time t = 500 Myr (final Stromgren sphere) for (left to right and top to bottom) C^-Ray, OTVET, CRASH, RSPH, ART, FTTE, SimpleX, FLASH-HC, 
and IFT. 



solving the non-equilibrium chemistry equations. These equations 
are stiff and thus generally require implicit solution methods. Such 
methods are generally expensive, however, so often certain approx- 
imations are used to speed-up the calculations. In order to test them, 
we performed the following simple test with a single, optically 
thin zone. We start with a completely neutral zone at time t = 0. 
We then applied photoionizing flux of F = 10^^ photons/s/cm^, 
with 10^ K black-body spectrum for 0.5 Myr, which results in the 
gas parcel becoming heated and highly ionized. Thereafter, the ion- 
izing flux is switched off and the zone cools down and recombines 
for further 5 Myr. The zone contains only hydrogen gas with num- 
ber density of n = 1 cm~^, and initial temperature of Ti — 100 K. 
Our results are shown in Figure |5] 

All codes agree very well in terms of the evolution of the neu- 
tral fraction (top panel), with the sole exception of SimpleX, in 
which case both the speed with which the gas parcel ionizes up and 
the achieved level of ionization are significantly different from the 
rest. The reason for this discrepancy is that currently this code does 
not solve the energy equation to find the gas temperature but has to 
assume a value instead (T = 10^ K in this test). FTTE finds slightly 
higher temperatures after its initial rise and correspondingly lower 
neutral fractions. 



Some differences are also seen after time t ^ 0.1 Myr, at 
which point there is a slight rise in temperature and corresponding 
dip in the neutral fraction. These occur around the time when the 
recombinations start becoming important, since tree ^ 0.1 Myr, 
which gives rise to slight additional heating. About half of the codes 
predict somewhat lower temperature rise then the rest. The cool- 
ing/recombination phase after source turn off demonstrates good 
agreement between the codes, although there is small difference in 
the final temperatures reached which is due to small differences in 
the hydrogen line cooling rates, resulting in slightly different tem- 
peratures at which the cooling becomes inefficient. 



3.2 Test 1: Pure-hydrogen isothermal H II region expansion 

This test is the classical problem of an H II re gion expansion in 
an uniform g as around a single ionizing source iStromgrenl Il939li 
Suitzer 1978). A steady, monochromatic Qiv — 13.6 eV) source 
emitting ionizing photons per unit time is turning on in an 
initially-neutral, uniform-density, static environment with hydro- 
gen number density uh. For this test we assume that the temper- 
ature is fixed at T = 10^ K. Under these conditions, and if we 
assume that the front is sharp (i.e. that it is infinitely-thin, with the 
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Figure 8. Test 1 (H II region expansion in an uniform gas at fixed temperature): Spherically-averaged profiles for ionized fractions x and neutral fractions 
xhi = 1 — X at times t = 30 Myr and 500 Myr vs. dimensionless radius (in units of the box size). 




Figure 7. Test 1 (H II region expansion in an uniform gas at fixed tempera- 
ture): The evolution of the position and velocity of the I-front. 



gas inside fully-ionized and the gas outside fully-neutral) there is a 
well-known analytical solution for the evolution of the I-front ra- 
dius, n, and velocity, vi, given by 

ri = rs [1 - exp(-t/trec)]^^^ , 
rs exp(-t/trec) 



vi 



where 



rs 



3tr, 



[l-exp(-t/trec)]'/' 



(5) 
(6) 



1 1/3 



47ras(T)n^ 



(7) 



is the Stromgren radius, i.e. the final, maximum size of the ionized 
region at which point recombinations inside it balance the incoming 
photons and the H II region expansion stops. The Stromgren radius 
is obtained from 



diUeriHOiBiT), 



(8) 



i.e. by balancing the number of recombinations with the number of 
ionizing photons arriving along a given line of sight (LOS). Here 
He is the electron density. 



[aB{T)nu] ^ , 



(9) 



is the recombination time, and a b (T) is the Case B recombination 
coefficient of hydrogen in the ionized region at temperature T. The 
H II region initially expands quickly and then slows considerably 
as the evolution time approaches the recombination time, t ^ tree, 
at which point the recombinations start balancing the ionizations 
and the H II region approaches its Stromgren radius. At a few re- 
combination times I-front stops at radius n = rs and in absence of 
gas motions remains static thereafter. The photon mean-free-path is 
given by 



Amfp — 



0.041 pc. 



(10) 



riHcro 

The particular numerical parameters we used for this test 
are as follows: computational box dimension L = 6.6 kpc, gas 
number density uh = 10 ~^ cm~^, initial ionization fraction 
(given by collisional equilibrium) x — 1.2 x 10~^, and ioniza- 
tion rate iV^ = 5 x 10^^ photons s~^. The source is at the cor- 
ner of the box). For these parameters the recombination time is 
tree = 3.86 X 10^^ s = 122.4 Myr. Assuming a recombination rate 
asiT) = 2.59 x 10"^^ cm^s-^ at T = 10^ K, then rs = 5.4 
kpc. The simulation time is tsim = 500 Myr ^ 4 tree The required 
outputs are the neutral fraction of hydrogen on the whole grid at 
times t = 10, 30, 100, 200, and 500 Myr, and the I-front position 
(defined by the 50% neutral fraction) and velocity vs. time along 
the X-axis. 
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Figure 9. Test 1 (H II region expansion in an uniform gas at fixed temperature): Fraction of cells with a given neutral fraction, xhi = 1 — x at times (left) 
t = 10 Myr, (middle) 100 Myr and (right) 500 Myr. 



In Figure |6l we show images of the neutral fraction in the 
z = plane at time t = 500 Myr, at which point the equilibrium 
Stromgren sphere is reached. The size of the final ionized region is 
in very good agreement between the codes. In most cases the H II 
region is nicely spherical, although some anisotropics exist in the 
CRASH and SimpleX results. In the first case these are due to the 
Monte-Carlo random sampling nature of this code, while in the sec- 
ond case it is due to the unstructured grid used by that code, which 
had to be interpolated on the regular grid format used for this com- 
parison. There are also certain differences in the H II region ionized 
structure, e.g. in the thickness of the ionized-neutral transition at the 
Stromgren sphere boundary. The inherent thickness of this transi- 
tion (defined as the radial distance between 0. 1 and 0.9 ionized frac- 
tion points) for monochromatic spectrum is « ISAmfp = 0.74 kpc, 
or about 14 simulation cells, equal to 11% of the simulation box 
size. This thickness is indicated in the upper left panel of Figure|6l 
Most codes find widths which are very close to this expected value. 
Only the OTVET, CRASH and SimpleX codes find thicker tran- 
sitions due to the inherently greater diffusivity of these methods, 
which spreads out the transition. For the same reason the highly- 
ionized proximity region of the source (blue-black colors) is no- 
ticeably smaller for the same two codes. 

In Figure we show the evolution of the I-front position and 
velocity. The analytical results in equation (|6j are shown as well 
(black, solid lines). All codes track the I-front correctly, with the 
position never varying by more than 5% from the analytical solu- 
tion. These small differences are partly due to differences in our 
recombination rates, as discussed above, and partly a consequence 
of our (somewhat arbitrary) definition of the I-front position as 
the point of 50% ionization. Our chosen parameters are such that 
the I-front internal structure is well-resolved, and the I-front in- 
trinsic thickness is larger than the discrepancies between the dif- 
ferent codes. The IFT code in particular tracks the I-front almost 
perfectly, as is expected for this code by construction. The ray- 
tracing codes agree between themselves a bit better than they do 
with the moment-based method OTVET. This is again related to 
the different, somewhat more diffusive nature of the last code. The 
I-front velocities also show excellent agreement with the analyt- 
ical result, at least until late times (at few recombination times), 
at which point the I-front essentially stops and its remaining slow 
motion forward is not possible to resolve with the relatively coarse 
resolution adopted for our test. The I-front at this point is moving 
so slowly that most of its remaining motion takes place within a 




Figure 10. Test 1 (H II region expansion in an uniform gas at fixed temper- 
ature): Evolution of the total neutral fraction. 



single grid cell for extended periods of time and thus falls below 
our resolution there. 

In Figure |8] we plot the spherically-averaged radial profiles of 
the ionized and the neutral fraction. In the left panel we show these 
profiles at t = 30 Myr, during the early, fast expansion of the I- 
front. Most of the ray-tracing codes (C^-Ray, ART, FLASH-HC 
and IFT) agree excellently at all radii. The OTVET, CRASH and 
Simplex codes appear more diffusive, finding a thicker I-front tran- 
sition and lower ionized fraction inside the H II region. The Zeus 
code also derives lower ionized fractions inside the H II region due 
to its slightly higher recombinational coefficient. The RSPH code is 
intermediate between the two groups of codes, finding essentially 
the same neutral gas profile inside the H II region as the ray-tracing 
codes, but a slightly thicker I-front, i.e. the ionized fraction drops 
more slowly ahead of the I-front. 
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Figure 11. Test 2 (H II region expansion in an uniform gas with varying temperature): Images of the H I fraction, cut through the simulation volume at 
coordinate ;2 = at time t = 10 Myr for (left to right and top to bottom) C^-Ray, OTVET, CRASH, RSPH, ART, FTTE, and IFT. 



The same differences persist in the ionized structure of the fi- 
nal Stromgren sphere at t = 500 Myr (Figure [Sj left panel). The 
majority of the ray-tracing codes again agree fairly well. The IFT 
code is based on the exact analytical solution of this particular prob- 
lem, and thus to a significant extent could be considered a substitute 
for the analytical H II region structure. Its differences from the ex- 
act solution are only close to the I-front, where the non-equilibrium 
effects dominate, while IFT currently assumes equilibrium chem- 
istry. Away from the I-front, however, the ionized state of the gas is 
in equilibrium and there all ray-tracing codes agree perfectly. The 
Simplex, OTVET and CRASH codes find thicker sphere bound- 
aries and lower ionized fractions inside, but the first two codes find 
slightly smaller Stromgren spheres, while the last finds a slightly 
larger one. The Zeus code finds lower ionized fractions inside the 
ionized region and a somewhat thicker I-front, but an overall H II 
region size that agrees with the other ray-tracing codes. The lower 
ionization is due to the current restriction of Zeus to monochro- 
matic radiative transfer, with its lower postfront temperatures and 
hence higher recombination rates. 

In Figure|9lwe show histograms of the fraction of cells with a 
given neutral fraction during the early, fast expansion phase (at time 
t = 10 Myr; left), when it starts slowing down (t = 100 Myr, close 
to one recombination time; middle), and when the final Stromgren 



sphere is reached (t = 500 Myr; right). These histograms reflect 
the differences in the I-front transition thickness and internal struc- 
ture. All codes predict a transitional region of similar size, which 
contains a few percent of the total volume. In detail, however, once 
again the results fall into two main groups. One group includes 
most of the ray-tracing codes, which agree perfectly at all times 
and predict thin I-fronts close to the analytical prediction. The 
other group includes the more diffusive schemes, namely OTVET, 
CRASH, RSPH and SimpleX, which all find somewhat thicker I- 
fronts. During the expansion phase of the H II region these three 
codes agree well between themselves, but they disagree somewhat 
on the structure of the final equilibrium Stromgren sphere, particu- 
larly in the proximity region of the source. The IFT code histograms 
differ significantly from the rest, due to its assumed equilibrium 
chemistry, which is not quite correct at the I-front. 



Finally, the evolution of the globally- averaged neutral frac- 
tions is shown in Fig ure[T0| The same trends are evident, with the 
ray-tracing codes agreeing closely among themselves, while the 
OTVET finds about 10% more neutral material at the final time, 
due to the different ionization structure obtained by this method. 
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Figure 12. Test 2 (H n region expansion in an uniform gas with varying temperature): Images of the temperature, cut through the simulation volume at 
coordinate ;2 = at time t = 10 Myr for (left to right and top to bottom) C^-Ray, OTVET, CRASH, RSPH, ART, FTTE, and IFT. 



3.3 Test 2: H II region expansion: the temperature state 

Test 2 solves essentially the same problem as Test 1, but the ion- 
izing source is assumed to have a 10^ K black-body spectrum and 
we allow the gas temperature to vary due to heating and cooling 
processes, as determined by the energy equation. The test geome- 
try and gas density are the same as in Test 1. The gas is initially 
fully neutral and has a temperature of T = 100 K. 

In Figure [n] we show images of the neutral fraction on the 
z = plane at time t = 10 Myr, during the initial fast expan- 
sion phase of the H II region. All of the results agree fairly well 
on the overall size of the ionized region and its internal structure. 
Again there are modest differences in the thickness of the I-front 
and the ionizing source proximity region, e.g. the CRASH code 
again produces a somewhat thicker transition and smaller proximity 
region. The IFT code finds significantly sharper I-fronts due to its 
equilibrium chemistry, but the internal H II region structure away 
from the front (which is close to equilibrium) and overall size of 
the ionized region both agree well with the rest of the codes. The 
temperature structure of the H II region, on the other hand, demon- 
strates significant differences (Figure[T2|. These stem largely from 
the different way the codes handle spectrum hardening, i.e. the long 
mean free paths of the high-energy photons due to the much lower 



photoionization cross-section at high frequencies. These long mean 
free paths result in a much thicker I-front transition and a signifi- 
cant pre-heating ahead of the actual I-front, since the high energy 
photons heat the gas, but there are not enough of them to ionize 
it. The CRASH code, which follows multiple bins in frequency, 
finds a larger pre-heated region than the other codes. There are also 
significant anisotropics in the CRASH results, due to the Monte- 
Carlo sampling method used (since not many high-energy photon 
packets are sent, leading to undersampling in angle). In production 
runs, a multifrequency treatment of the single photon packets has 
been introduced, reducing the anisotropics in the results. The tem- 
perature results of C^-Ray, RSPH and ART codes agree fairly well 
among themselves, while OTVET and FTTE give much less spec- 
trum hardening. Finally, IFT assumes the I-front is sharp, and does 
not have spectrum hardening by construction. 

The same trends persist at later times, when the I-front is ap- 
proaching the Stromgren sphere (Figures fT3l andfT4l. Once again 
the H II regions predicted by all the codes are similar in size and in- 
ternal structure, but with a little different I-front thickness in terms 
of neutral fraction and significant differences in terms of spectral 
hardening. The FTTE still gives a very sharp I-front, while OTVET 
finds somewhat less hardening, but its later-time result is more sim- 
ilar to the other codes than at early times. 
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Figure 13. Test 2 (H II re^ 
coordinate z = at time t - 



;ion expansion in an uniform gas with varying temperature): Images of the H I fraction, cut through the simulation volume at 
= 100 Myr for (left to right and top to bottom) C^-Ray, OTVET, CRASH, RSPH, ART, FTTE, and IFT. 



In Figure[^we plot the the position and velocity of the I-front 
vs. time. Unlike Test 1, in this case there is no closed-form ana- 
lytical solution since the recombination coefficients vary with the 
spatially-varying temperature. Nevertheless, as a point of reference 
we have again shown the analytical solution in equations (as- 
suming T — 10^ K). All the codes find slightly larger H II regions 
and slightly faster I-front propagation compared to this analytical 
solution. This is to be expected due to the temperature being higher 
than 10^ K and the inverse temperature dependence of the recombi- 
nation coefficient. The C^-Ray, RSPH and FTTE results agree per- 
fectly among themselves, to 1%, as do the results from OTVET, 
ART, and Zeus, again among themselves. These two groups of re- 
sults differ by 10%, however, while CRASH and IFT find an 
H II region size intermediate between the two groups. 

In Figure[l6lwe show the spherically-averaged radial profiles 
of the neutral and ionized fractions during the fast expansion phase 
(t = 10 Myr, left), the slowing-down phase (t = 100 Myr, middle) 
and the final Stromgren sphere (t = 500 Myr, right). These con- 
firm, in a more quantitative way, the trends already noted based on 
the 2D images above. The profiles from C^-Ray and RSPH codes 
are in excellent agreement at all radii and all times. The IFT code 
closely agrees with them in the source proximity region, where the 
gas ionized state is at equilibrium, but diverges around the I-front 



(due to its assumed equilibrium chemistry) and ahead of the I-front 
(due to its assumption that the front is sharp). Compared to these 
codes, CRASH and ART find slightly higher neutral fractions close 
to the ionizing source, but agree well with C^-Ray and RSPH ahead 
of the I-front in the spectral hardening region. The FTTE code is 
in excellent agreement with C^-Ray and RSPH codes close to the 
source, but its I-front is much sharper. The OTVET code also finds 
a somewhat sharper I-front, but to a much lesser extent than the 
FTTE code, while close to the source its neutral fraction is only 
slightly higher than the majority of codes, and in close agreement 
with the ART code. Finally, the I-front derived by the Zeus code is 
very sharp. This is due to the use of only single-energy photons by 
this code, which does not allow for spectral hardening. 

The corresponding spherically-averaged radial temperature 
profiles at the same three I-front evolutionary phases are shown 
in Figure [H] All results (except the one from Zeus, due to the 
monochromatic spectrum it used) agree well inside the ionized re- 
gion, with the differences arising largely due to slight differences 
in the cooling rates adopted. The more diffusive OTVET code does 
not show as sharp temperature rise in the source proximity as the 
other codes, but elsewhere the temperature structure it finds agrees 
with the majority of the codes. Again at the I-front and ahead of 
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it the differences between the results are significant, reflecting the 
different handUng of hard photons by the codes. 

In Figure fTsI we show the histograms of the fraction of cells 
with a given ionized fraction x at the same times as the radial pro- 
files above. During the early, fast expansion phase all codes agree 
well except the OTVET code, which finds a slightly thinner I-front, 
but an otherwise same histogram distribution shape, and CRASH, 
whose I-front is a bit thicker. IFT finds a different ionized fraction 
distribution, again as a consequence of its equilibrium chemistry, 
which is not correct at the I-front transition. Later, when the I-front 
slows down (t = 100 Myr) the same trends hold, but in addition 
the FTTE results start diverging significantly from the rest, finding 
notably smaller ionized region and a quite different shape distri- 
bution in the largely-neutral regions. This reflects its much sharper 
I-front with little spectrum hardening, as noted above. Finally, the 
ionized fraction histogram corresponding to the Stromgren sphere 
(T = 500 Myr) shows similar differences. The C^-Ray, ART and 
RSPH codes again agree very closely, and CRASH also finds a 
similar distribution, but with a thicker I-front and correspondingly 
fewer neutral cells. The OTVET distribution follows a roughly sim- 
ilar shape but with a thinner front transition and more neutral cells, 
while FTTE agrees well with the other ray-tracing codes in the 



highly-ionized region, but still diverges considerably at the I-front 
and ahead of it. 

The corresponding histograms of temperature are shown in 
Figure [191 The ionized gas temperatures found by all codes have 
a strong peak at slightly above 10^ K, with only small variations 
in the peak position between the different codes. This peak was to 
be expected, as a consequence of the combination of photoheating 
and hydrogen line cooling (which peaks around 10"^ K). At higher 
temperatures (corresponding to cells close to the ionizing source) 
some differences emerge. The ray-tracing codes (C^-Ray, CRASH, 
ART, RSPH, IFFT, IFT) largely agree among themselves, with only 
C^-ray finding slightly larger fraction of hot cells. OTVET, on 
the other hand does not predict any cells with temperature above 
^ 16, 000 K, due to missing the very hot proximity region of the 
source, as noted above. Below ~ 8, 000 K, on the other hand, the 
differences between results are more significant, reflecting the vari- 
ations in the I-front thickness and spectral hardening noted above. 
The results from CRASH, ART and RSPH agree well at all times, 
while the C^-Ray histograms have similar shape, but with some 
offset, due to its current somewhat simplified handling of the en- 
ergy input which uses a single bin in frequency. The OTVET, FTTE 
and IFT codes find much smaller pre-heating ahead of the I-front, 
and thus different distributions. 
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Figure 15. Test 2 (H II region expansion in an uniform gas with varying temperature): The evolution of the position and velocity of the I-front. 
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Figure 16. Test 2 (H II region expansion in an uniform gas with varying temperature): Spherically-averaged ionized fraction x and neutral fraction 1 — x 
profiles at times t = 10 Myr, 100 Myr and 500 Myr. 
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Figure 17. Test 2 (H II region expansion in an uniform gas with varying temperature): Spherically-averaged temperature profiles at times t = 10 Myr, 100 
Myr and 500 Myr. 
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Figure 18. Test 2 (H II region expansion in an uniform gas with varying temperature): Fraction of cells with a given ionized fraction, x, at times (left) t = 10 
Myr, (middle) 100 Myr and (right) 500 Myr. 
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Figure 19. Test 2 (H II region expansion in an uniform gas with varying temperature): Fraction of cells with a given temperature T at times (left) t = 10 Myr, 
(middle) 100 Myr and (right) 500 Myr. 
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Figure 21. Test 3 (I-front trapping in dense clump): The evolution of the position and velocity of the I-front. 
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Figure 20. Test 2 (H II region expansion in an uniform gas with varying 
temperature): Evolution of the total neutral gas fraction. 



Finally, in Figure |20| we show the evolution of the total neu- 
tral gas fraction. All codes agree well on the final neutral fraction, 
within 25% or better. The differences are readily understood in 
terms of the different recombination rates, mostly as a consequence 
of the somewhat different temperatures found inside the H II re- 
gion, in addition to the small differences in the recombination rate 
fits used. 



3.4 Test 3: I-front trapping in a dense dump and the 
formation of a shadow. 

Test 3 examines the propagation of a plane-parallel I-front and its 
trapping by a dense, uniform, spherical clump. The condition for 
an I-front to be trapped by a clump of gas with n umber density uh 
can be derived as follows ('Shapiro etalll2004 . Let us define the 
Stromgren length Is (r) at impact parameter r from the clump cen- 
ter using equation Js}, but in this case following lines-of- sight for 
each impact parameter. We then can define the "Stromgren num- 
ber" for the clump as Ls = 2rciump/^s(0), where rdump is the 
clump radius and is(0) is the Stromgren length at zero impact pa- 
rameter. Then, if > 1 the clump is able to trap the I-front, while 
if Ls < 1, the clump would be unable to trap the I-front and instead 
would be flash-ionized by its passage. 

For a uniform-density clump equation (|8} reduces to 



(11) 



(2) 2 ■ 



and the Stromgren number is given by 

J 2rciumpas(T)n^ 
Ls — ^ . 



(12) 



The numerical parameters for Test 3 are as follows: the spec- 
trum is a black-body with effective temperature Te/ / = 10^ K and 
constant ionizing photon flux, F — 10^ s~^cm~^, incident to the 
2/ = box side; the hydrogen number density and initial tempera- 
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Figure 22. Test 3 (I-front trapping in dense clump): Images of the H I fraction, cut through the simulation volume at midplane at time t = 1 Myr for C^-Ray, 
CRASH, RSPH, FTTE, Flash-HC, IFT and Coral. 



ture of the environment are riout = 2 x 10~^ cm~^ and Tout,init = 
8, 000 K, while inside the clump they are riciump = 200nout = 
0.04 cm~^ and Tciump,init = 40 K. The box size is L = 6.6 
kpc, the radius of the clump is rdump = 0.8 kpc, and its center 
is at (xc, Vc, Zc) = (5, 3.3, 3.3) kpc, or (xc, yc, Zc) = (97, 64, 64) 
cells, and the evolution time is 15 Myr. For these parameters and 
assuming for simplicity that the Case B recombination coefficient 
is given by as(T) = 2.59 x 10"^^ (T/10^ ii^)"^/^, we obtain 
Is ^ 0.78(T/10^ i^)^/^ kpc, and ^ 2.05(T/10^ i^)"^/^; 
thus, along the axis of symmetry the I-front should be trapped ap- 
proximately at the center of the clump for T = 10^ K. In reality, 
the temperature could be expected to be somewhat different and 
spatially-varying, but to a rough first approximation this estimate 
should hold. 

In fact, the I-front does get trapped as expected, slightly be- 
yond the clump center. In Figure|2]we plot (n — Xc)/is, the evo- 
lution of the position of the I-front with respect to the clump center 
in units of the Stromgren length (top panel) and the corresponding 
velocity evolution (in units of 2c5,i(T = lO^K) = 2(p/p)^/^ 
twice the isothermal sound speed in gas at temperature of 10^ K), 
both vs. t/trec,o, time in units of the recombination time inside the 
clump (which is ^ 3 Myr at 10^ K). The I-front is initially highly 
supersonic due to the low density outside the clump. Once it enters 
the clump it shows down sharply, to about 20 times the sound speed, 
by the same factor as the density jump at the clump boundary. As 



it penetrates further into the clump it approaches its (inverse, i.e. 
outside-in) Stromgren radius at time t ^ trec,o, at which point the 
propagation slows down even further until the I-front is trapped 
after a few recombination times. The velocity drops below 2cs,/, 
at which point if gas motions were allowed the I-front would be- 
come slow D-type (i.e. coupled to the gas motion, rather than much 
faster than them). All codes capture these basic phases of the trap- 
ping process correctly and agree well on both the front position and 
velocity. We note here that the FLASH-HC code currently does not 
have the ability to track fast I-fronts, so its data starts only after the 
front has slowed down. The IFT method assumes a sharp front, and 
thus does not allow pre-heating and partial ionization ahead of the 
front, which results in its being slowed down more abruptly than 
is the case for the other results. Due to some diffusion, the RSPH 
code finds that the front slows down slightly before the I-front actu- 
ally enters the clump. There are also minor differences in the later 
stages of the evolution, to be discussed in more detail below. 

In Figure |22l we show the images of the neutral gas fraction 
on the plane through the centre of the clump at time t = 1 Myr, 
when the I-front is already inside the clump, but still not trapped 
and moving supersonically. The ionizing source is far to the left 
of the box. All results show a sharp shadow behind the clump, as 
expected for such a dense, optically-thick clump. Only the RSPH 
code shows diffusion at the shadow boundaries, due to the intrin- 
sic difficulty of representing such a sharply-discontinuous density 
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Figure 23. Test 3 (I-front trapping in dense clump): Images of the temperature, cut through the simulation volume at midplane at time t = 1 Myr for C^-Ray, 
CRASH, RSPH, FTTE, Flash-HC, IFT and Coral. 



distribution with SPH particles and the corresponding smoothing 
kernel. The FLASH-HC code derives a noticeably sharper I-front 
in both the clump and the external medium (where it is the only 
result to still have some gas with neutral fraction above ^ 10~^). 
This is due to its current inability to correctly track fast I-fronts, 
as discussed above, which leads to somewhat incorrect early evolu- 
tion. The corresponding temperature image cuts (Figure l23l show 
the same trends: FLASH-HC and IFT find a very sharp transition, 
while the rest of the codes agree reasonably well, with only minor 
differences in the pre-heating region. 

In Figure |24l we show the images of the neutral fraction at the 
final time of the simulation, t = 15 Myr. All codes except CRASH 
find very similar ionized structure inside the clump. The CRASH 
result has significantly higher neutral fraction inside the clump, and 
correspondingly larger shadow behind the clump, as well as slightly 
higher neutral gas fraction in the low-density gas. This could be 
due to the fact that, as mentioned in § 13.31 CRASH follows multi- 
ple bins in frequency over a wider frequency range with respect to 
the other codes; this results in a higher ionizing power at high fre- 
quencies, which also have smaller photo-ionization cross-sections. 
This in turn could be the origin of the lower ionization state of the 
clump and of the low-density gas. The RSPH result again exhibits 
significant diffusion around the edges of the shadow. The corre- 
sponding temperature structures, on the other hand, (Figure l25l 
show some differences, which stem from the different treatments 



of the energy equation and spectral hardening by the codes. The 
FTTE and IFT codes find almost no pre-heating in the shielded re- 
gion and the shadow behind it. C^-Ray and CRASH get smaller 
self- shielded regions and some hard photons penetrating into the 
sides of the shadow. Finally, RSPH, FLASH-HC and Coral find al- 
most no gas that is completely self- shielded, but still find sufficient 
column densities to create temperature- stratified shadows similar to 
the ones found by C^-Ray and CRASH codes, albeit at higher tem- 
perature levels. The Coral result also has a thin, highly-heated shell 
at the source side of the clump, resulting from this code's prob- 
lems in properly finding the temperature state in the first dense, 
optically-thick cells encountered by its rays, which leads to their 
overheating. In production runs this problem was corrected by in- 
creasing the resolution and decreasing the cell size so that cells are 
not as optically-thick, and by the gas motions, which quickly cool 
the gas down as it expands. The low-density gas outside the clump 
is somewhat cooler in the CRASH and RSPH results compared to 
the other codes. 

These observations are confirmed by the evolution of the mean 
ionized fraction and the mean temperature inside the clump shown 
in Figure l26l All codes agree very well on the evolution of the mean 
ionized fraction, except for CRASH, which finds about ^ 25% 
lower final ionized fraction, and for FLASH-HC, which early-on 
finds a lower ionized fraction, but catches up with the majority of 
the codes as the I-front becomes trapped. In terms of mean temper- 
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Figure 24. Test 3 (I-front trapping in dense clump): Images of the H I fraction, cut through the simulation volume at midplane at time f = 15 Myr for C^-Ray, 
CRASH, RSPH, FTTE, Flash-HC, IFT and Coral. 



ature, Coral, and to a lesser extent FLASH-HC find higher mean 
temperature due to the overheating of some cells mentioned above. 

In Figure]^ we show the ionized and neutral fraction pro- 
files along the axis of symmetry at three stages of the evolution - 
early {t — 1 Myr), during the slow-down due to recombinations 
(t = 3 Myr, about one recombination time in the clump) and late 
(t = 15 Myr). Only the region inside and around the clump is 
plotted in order to show details. In the pre-ionization, spectrum- 
hardening zone ahead of the main I-front all profiles agree fairly 
well at all times, except that the IFT profiles have a sharp I-front 
and no hardening by definition, and the FTTE current method ap- 
pears to produce no hardening, either. Otherwise, these two codes 
agree well with the others in the post-front region. Some differ- 
ences emerge in the position of the I-front (defined as the point of 
50% ionized fraction) and the neutral fraction profiles behind the 
I-front. During the initial, fast propagation phase of the front all 
codes agree on its position. However, CRASH and, to a lesser ex- 
tent, FLASH-HC both find consistently higher neutral fractions in 
the ionized part of the clump than the rest of the codes. As a conse- 
quence, they show that the I-front is trapped closer to the surface of 
the dense clump. The reason for this can be seen in the correspond- 
ing temperature profiles (Figure l29l. Both CRASH and FLASH- 
HC obtain a slightly lower temperature for the dense ionized gas, 
resulting in a higher recombination rate there. The RSPH and IFT 
results show some diffusion at the source- side clump boundary (at 



r/Lho^ ^ 0.64), resulting in a less sharp transition there, regard- 
less of the sharp discontinuity of the gas density. The temperatures 
in the post-front fraction of the clump otherwise agree quite well, 
with the exception of the first cells on the source side, where Coral 
finds very high temperatures, as was already mentioned. In the pre- 
front region the temperature profile results also agree fairly well. 
Only IFT and FTTE differ there, again producing a very sharp I- 
front. C^-Ray finds a bit less pre-heating in this region due to its 
single frequency bin method. Until t = 3 Myr all results except 
the FLASH-HC one agree that the shadow region right behind the 
clump is completely shielded and remains at the initial temperature 
(8,000 K). At the final time, however, more differences emerge. 
At this time only IFT and FTTE still find no pre-heating, and C^- 
Ray and CRASH find only little heating of the shadow. The Coral, 
RSPH and FLASH-HC results all find significant pre-heating, al- 
beit at different levels. RSPH is the only code that at the final time 
finds no partially shielded gas at the back of the clump, the presence 
of which is indicated by the temperature dip at r/Lbox 0.88 for 
all the other results. 

The final results we show for this test are histograms of the 
fraction of cells inside the clump with a given ionized fraction (Fig- 
ure |29j and with a given temperature (Figure l30t at the same three 
evolutionary stages discussed above. These histograms reflect the 
differences in the thickness and the internal structure of the I-fronts 
inside the clump. During the initial fast propagation and the slow- 
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Figure 25. Test 3 (I-front trapping in dense clump): Images of the temperature, cut through the simulation volume at midplane at time t = 15 Myr for C^-Ray, 
CRASH, RSPH, FTTE, Flash-HC, IFT and Coral. 



down phases the RSPH and Coral codes consistently find somewhat 
thicker I-front transitions than the other codes. The results of C^- 
Ray, CRASH and FLASH-HC start with the front somewhat thin- 
ner than the ones for RSPH and Coral, but the distribution changes 
as the I-front gets trapped. At the final time, the ionized fraction 
distributions are very similar, with only the one for CRASH be- 
ing slightly thicker. Finally, FTTE and IFT find significantly thin- 
ner front transitions and significant self-shielded gas fractions, as 
was noted before. The corresponding temperature distributions all 
show a strong peak at a similar temperature, a few tens of thou- 
sands of degrees, typical for gas heated by photoionization. This 
temperature is well above 10^ K because of the hot black-body 
spectrum of the source. At temperatures lower than this peak value, 
which largely correspond to the pre-heated zone ahead of the I- 
front, there again is broad agreement, apart from the FTTE and 
IFT codes which have a sharp front and little pre-heating. At the 
high-temperature end of the distribution C^-Ray, FTTE, FLASH- 
HC and Coral diverge from the rest, by finding a small fraction of 
very hot cells, while the majority of codes find essentially no cells 
hotter than the distribution peak. It should be noted, however, that 
the fraction of these cells is only 0.1-1% of the total. 



3.5 Test 4: Multiple sources in a cosmological density field 

Test 4 involves the propagation of I-fronts from multiple sources 
in a static cosmological density field. The initial condition is pro- 
vided by a time- slice (at redshift z = 9) from a cosmological N- 
body and gas-dynamic simulatio n performed u sing the cosmolog- 
ical PM-hTVD code by D. Ryu iRvu etalJ[T99 31 The simulation 
box size is 0.5 comoving Mpc, the resolution is 128^ cells, 
2 X 64^ particles. The halos in the simulation box were found using 
friends-of-friends halo finder with linking length of 0.25. For sim- 
plicity the initial temperature is fixed at T = 100 K everywhere. 
The ionizing sources are chosen so as to correspond to the 16 most 
massive halos in the box. We assume that these have a black-body 
spectrum with effective temperature Tgff = 10^ K. The ionizing 
photon production rate for each source is constant and is assigned 
assuming that each source lives ts = 3 Myr and emits = 250 
ionizing photons per atom during its lifetime, hence 



Qorrints 



(13) 



where M is the total halo mass, = 0.27, = 0.043, h = 0.7. 

For simplicity all sources are assumed to switch on at the same 
time. The boundary conditions are transmissive (i.e. photons leav- 
ing the computational box are lost, rather than coming back in as in 
periodic boundary conditions). The evolution time is t = 0.4 Myr. 
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Figure 26. Test 3 (I-front trapping in dense clump): The evolution of the ionized fraction (top) and the mean temperature (bottom) inside the dense clump. 




Figure 27. Test 3 (I-front trapping in dense clump): Line cuts of the ionized and neutral fraction along the axis of symmetry through the center of the clump 
at times t = 1 Myr, 3 Myr and 15 Myr. 
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Figure 28. Test 3 (I-front trapping in dense clump): Line cuts of the temperature along the axis of symmetry through the center of the clump at times t = 1 
Myr, 3 Myr and 15 Myr. 
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Figure 29. Test 3 (I-front trapping in dense clump): Histograms of the ionized fraction inside the clump at times t = 1,3 and 15 Myrs. 



In Figure lsTl we show slices of the H I fraction cut through the 
simulation box at coordinate z = ^box/2 and time t = 0.05 Myr, 
and in Figure |32l we show the corresponding temperature distri- 
butions the same time. SimpleX does not appear in the tempera- 
ture maps since currently the code does not follow the temperature 
evolution self-consistently, but instead assumes a constant temper- 
ature. In Figures |33l and |34l we show the same as above but at time 
t = 0.2 Myr. Some discrepancies are already evident from a visual 
inspection, which shows somewhat different morphologies, but still 
general agreement. In Figure |35l (left panel) we present the tem- 
poral evolution of the volume- (thin lines) and mass- (thick lines) 
weighted ionized fractions. While CRASH and FTTE find com- 
parable ionized fractions, C^-Ray and SimpleX produce slightly 
higher and lower values, respectively The lower value in SimpleX 
is obtained as a consequence of the temperature being fixed to 
10^ K, which is a little lower than that obtained by the other codes, 
resulting in a higher recombination rate. 

The finer sampling in CRASH of the high energy tail of the 
spectrum allows a higher resolution of its hardening. As the total 
energy is distributed differently from the other codes (more energy 



is in the hard photons), this results in less ionization/heating closer 
to the sources and more ionization/heating further away and a lower 
mean ionized fraction. 

In Figuresl36landl37lwe show histograms of the ionized frac- 
tion and the temperature at times t = 0.05, 0.2 and 0.4 Myr are 
shown. While C^-Ray and FTTE agree quite well, especially at 
later times, CRASH, as explained above, produces a thicker I-front 
transitions. The thickness of the I-front as found by SimpleX ap- 
pears to oscillate with time, initially starting thick, then becoming 
thinner, and thick again towards the end of the simulation. This 
might have been caused by the interpolation of its unstructured grid 
to the regular grid required for the results. In terms of temperature, 
CRASH finds a systematically lower peak 1.6 x 10^ K com- 
pared to - 6.3 X 10^ K for C^-Ray and FTTE, which agree well 
there) inside the ionized regions and a higher value in the lower 
density regions. This is once again due to the spectrum hardening 
effects discussed above. The peak produced by FTTE at very low 
temperatures arises because of the lack of spectrum hardening and 
sharp I-fronts consistently produced by this code, in which case no 
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Figure 30. Test 3 (I-front trapping in dense clump): Histograms of the temperature inside the clump at times t = !,?> and 15 Myrs. 



photons propagate ahead of the I-front and thus the gas away from 
the sources remains cold and neutral. 

To understand whether the differences discussed above can be 
in fact attributed to photon hardening and the different treatment 
of the high energy tail of the spectrum, we have repeated the same 
simulations with a softer black-body spectrum with effective tem- 
perature Teff = 3 X 10^ K, in which case the spectrum harden- 
ing and I-front spreading should be minimized. We note that this 
part was not done by the SimpleX code since currently it does not 
treat different spectra. Our results are shown in Figures |35l (right 
panel), |38]and|39l The averaged ionized fraction produced by C^- 
Ray and CRASH now have a much better agreement, as far fewer 
hard photons are present. However, FTTE still produces a lower 
value, although by only ^ 5%. CRASH still obtains a somewhat 
thicker ionizing front, due to the inherent adopted method, but the 
agreement is now better, especially at later times, as high energy 
photons are not present in this case. This is even more evident from 
an analysis of the temperature behaviour, where the agreement be- 
tween C^-Ray and CRASH is now very good at all times, while 
FTTE consistently finds higher temperatures inside the ionized re- 



gions. These higher temperatures seem to self-contradict the lower 
ionized fractions found by this code, while the recombination rate 
used by this code is consistent with the others, indicating a possible 
(modest) problem with photon conservation. 



4 SUMMARY AND CONCLUSIONS 

We have presented a detailed comparison of a large set of cosmo- 
logical radiative transfer methods on several common tests. The 
participating codes represent the full variety of existing methods, 
multiple ray-tracing and one moment method, which solve the ra- 
diative transfer on regular, adaptive or unstructured grids, even with 
no grid at all, but instead using particles to represent the density 
field. The comparison is a collaborative project involving most of 
the cosmological RT community. The results from this comparison 
will be publicly available for testing of future codes during their 
development. 

We began by comparing our basic physics, like chemistry, 
cooling rates and photoionization cross-sections, which came from 
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Figure 31. Test 4 (reionization of a cosmological density field): Images of the H I fraction, cut through the simulation volume at coordinate z = 2:box/2 
and time t = 0.05 Myr for C^-Ray (left), CRASH (middle), FTTE (right), and SimpleX (bottom). The black-body spectrum has an effective temperature 
Teff = 105 K. 






Figure 32. Test 4 (reionization of a cosmological density field): Images of the temperature, cut through the simulation volume at coordinate z = 2;box/2 and 
time t = 0.05 Myr for C^-Ray (left), CRASH (middle), and FTTE (right). The black-body spectrum has an effective temperature Teff = 10^ K. 



a variety of sources, and evaluated the effects these have on the 
propagation of I-fronts. We concluded that even quite approximate 
rates generally result in relatively modest divergences in the results. 
The discrepancies were limited to ~ 5 — 6% in I-front position 
but were somewhat larger in velocity and the internal structure of 
the H II region (but never exceeding 20-40%, and usually much 
smaller). 

Then, we turned to some simple, but instructive and 
cosmologically-interesting problems which tested the different as- 
pects encountered in realistic applications. In summary, the results 
showed that at fixed temperature and for monochromatic ionizing 
spectrum all of these methods track I-fronts well, to within a few 
per cent accuracy. There are some differences in the inherent thick- 
ness (due to finite mean free path) and the internal structure of the 
I-fronts, related to some intrinsic diffusivity of some of the meth- 
ods. 

Somewhat greater differences emerge when the temperature 
is allowed to vary, although even in this case the agreement is very 
good, typically within 10-20% in the ionized fraction. Other varia- 
tions between the results are due to the variety of multi-frequency 
radiative transfer methods and consequent spectrum hardening and 



pre-heating ahead of the I-fronts. Some of the algorithms consis- 
tently find quite sharp I-fronts, with little spectrum hardening and 
pre-heating, while the majority of the codes handle these features 
more precisely and in reasonable agreement among themselves. 



We conclude that the various approximations employed by the 
different methods perform quite well and generally produce con- 
sistent and reliable results. There are, however, certain differences 
between the available methods and the most appropriate method 
for any particular problem would vary. Therefore, the code em- 
ployed should be chosen with care depending on the specific ques- 
tions whose answer is sought. For example, for problems in which 
the I-front structure is not resolved (due to coarse simulation res- 
olution compared to the characteristic front width), most methods 
considered here should perform well and the main criterion should 
be their numerical efficiency. On the other hand, in problems where 
the I-front width and spectral hardening matters significantly, atten- 
tion should be paid to the accurate multi-frequency treatment and 
the energy equation. 
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Figure 33. Test 4 (reionization of a cosmological density field): Images of the H I fraction, cut through the simulation volume at coordinate z = z\yo^/2 
and time t = 0.2 Myr for C^-Ray (left), CRASH (middle), FTTE (right), and SimpleX (bottom). The black-body spectrum has an effective temperature 
Teff = 105 K. 




Figure 34. Test 4 (reionization of a cosmological density field): Images of the temperature, cut through the simulation volume at coordinate z = z\yo 
time t = 0.2 Myr for C^-Ray (left), CRASH (middle), and FTTE (right). The black-body spectrum has an effective temperature Teff = 10^ K. 
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Figure 35. Test 4 (reionization of a cosmological density field): The evolution of the volume- and mass-weighted ionized fractions, Xv (thin lines) and Xm 
(thick lines), for a black-body source spectra with Teff = 10^ K (left) and Tgff = 3 x 10^ K (right). 
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Figure 36. Test 4 (reionization of a cosmological density field): Histograms of the neutral fraction at times t = 0.05, 0.2 and 0.4 Myrs for a black-body source 
spectrum with Teff = 10^ K. 
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Figure 37. Test 4 (reionization of a cosmological density field): Histograms of the temperature at times t = 0.05, 0.2 and 0.4 Myrs for a black-body source 
spectrum with Tgff = 10^ K. 
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Figure 38. Test 4 (reionization of a cosmological density field): Histograms of the ionized fraction at times t = 0.05, 0.2 and 0.4 Myrs for a black-body 
source spectrum with Tgff = 3 x 10^ K. 
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Figure 39. Test 4 (reionization of a cosmological density field): Histograms of the temperature at times t = 0.05, 0.2 and 0.4 Myrs for a black-body source 
spectrum with Tgff = 3 x lO'^ K. 
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